
!------------------------------------------------------------------------!
!  The Community Multiscale Air Quality (CMAQ) system software is in     !
!  continuous development by various groups and is based on information  !
!  from these groups: Federal Government employees, contractors working  !
!  within a United States Government contract, and non-Federal sources   !
!  including research institutions.  These groups give the Government    !
!  permission to use, prepare derivative works of, and distribute copies !
!  of their work in the CMAQ system to the public and to permit others   !
!  to do so.  The United States Environmental Protection Agency          !
!  therefore grants similar permission to use the CMAQ system software,  !
!  but users are requested to provide copies of derivative works or      !
!  products designed to operate in the CMAQ system to the United States  !
!  Government without restrictions as to use by others.  Software        !
!  that is used with the CMAQ system but distributed under the GNU       !
!  General Public License or the GNU Lesser General Public License is    !
!  subject to their copyright restrictions.                              !
!------------------------------------------------------------------------!

       SUBROUTINE HRG3( DTC )

C**********************************************************************
C
C  FUNCTION:  To solve for the concentration of C2O3 and PAN
C
C  PRECONDITIONS: For the CRACMM1AMORE_AQ mechanism
C
C  KEY SUBROUTINES/FUNCTIONS CALLED: None
C
C  REVISION HISTORY: Created by EBI solver program, Jun 16, 2022
C
C   18 Jul 14 B.Hutzell: revised to use real(8) variables
C**********************************************************************
      USE HRDATA

      IMPLICIT NONE

C..INCLUDES:  NONE

C..ARGUMENTS:
      REAL( 8 ), INTENT( IN ) :: DTC              ! Time step


C..PARAMETERS: NONE


C..EXTERNAL FUNCTIONS: NONE


C..SAVED LOCAL VARIABLES:
!     CHARACTER( 16 ), SAVE  :: PNAME = 'HRG3'      ! Program name


C..SCRATCH LOCAL VARIABLES:
      REAL( 8 ) ::   A, B, C, Q   ! Quadratic equation terms
      REAL( 8 ) ::   CMN          ! Temp scalar
      REAL( 8 ) ::   L8           ! Loss of CCO_O2
      REAL( 8 ) ::   L9           ! Loss of PAN
      REAL( 8 ) ::   P8           ! Production of CCO_O2

      REAL( 8 ) ::   K8_8         ! Kaco3+aco3 * delta t
      REAL( 8 ) ::   R8_9         ! Kpan-->aco3 * delta t
      REAL( 8 ) ::   R9_8         ! Kaco3+no2-->pan * [NO2] * delta t

C**********************************************************************


c..Production of ACO3 (except from PAN )
      P8 =                 RXRAT(    14 )      ! ACT=ACO3+MO2
     &   +    7.8400D-01 * RXRAT(    16 )      ! UALD=0.7840D+00*ACO3+...
     &   +    9.0000D-01 * RXRAT(    19 )      ! MEK=0.9000D+00*ACO3+ETHP+...
     &   +    5.0000D-01 * RXRAT(    20 )      ! KET=0.5000D+00*ACO3+...
     &   +                 RXRAT(    21 )      ! HKET=ACO3+HO2+HCHO
     &   +    6.7000D-01 * RXRAT(    22 )      ! MACR=0.6700D+00*ACO3+...
     &   +                 RXRAT(    27 )      ! MGLY=ACO3+HO2+CO
     &   +    2.5000D-01 * RXRAT(    28 )      ! DCB1=0.2500D+00*ACO3+...
     &   +    2.5000D-01 * RXRAT(    29 )      ! DCB2=0.2500D+00*ACO3+...
     &   +    1.0000D-01 * RXRAT(    88 )      ! ISO+O3=0.1000D+00*ACO3+...
     &   +    6.0000D-02 * RXRAT(    99 )      ! ISHP=0.6000D-01*ACO3+...
     &   +                 RXRAT(   109 )      ! ACD+HO=ACO3
     &   +    3.1300D-01 * RXRAT(   117 )      ! UALD+HO=0.3130D+00*ACO3+...
     &   +                 RXRAT(   119 )      ! MGLY+HO=ACO3+CO
     &   +    6.5000D-01 * RXRAT(   138 )      ! PAA+HO=0.6500D+00*ACO3+...
     &   +    9.0000D-02 * RXRAT(   146 )      ! OLI+O3=0.9000D-01*ACO3+...
     &   +    1.0000D-01 * RXRAT(   151 )      ! MACR+O3=0.1000D+00*ACO3+...
     &   +    2.8000D-01 * RXRAT(   152 )      ! MVK+O3=0.2800D+00*ACO3+...
     &   +    2.0000D-03 * RXRAT(   153 )      ! UALD+O3=0.2000D-02*ACO3+...
     &   +                 RXRAT(   165 )      ! ACD+NO3=ACO3+HNO3
     &   +                 RXRAT(   170 )      ! MGLY+NO3=ACO3+CO+HNO3
     &   +                 RXRAT(   206 )      ! ACTP+NO=ACO3+NO2+HCHO
     &   +    2.3000D-01 * RXRAT(   208 )      ! KETP+NO=0.2300D+00*ACO3+...
     &   +    3.5000D-01 * RXRAT(   209 )      ! MACP+NO=0.3500D+00*ACO3+...
     &   +    7.0000D-01 * RXRAT(   211 )      ! MVKP+NO=0.7000D+00*ACO3+...
     &   +    1.5000D-01 * RXRAT(   248 )      ! ACTP+HO2=0.1500D+00*ACO3+...
     &   +    5.0000D-01 * RXRAT(   285 )      ! ACTP+MO2=0.5000D+00*ACO3+...
     &   +    2.6900D-01 * RXRAT(   288 )      ! MACP+MO2=0.2690D+00*ACO3+...
     &   +    1.1600D+00 * RXRAT(   290 )      ! MVKP+MO2=0.1160D+01*ACO3+...
     &   +    1.6000D-01 * RXRAT(   327 )      ! MVKP+ACO3=0.1160D+01*ACO3+...
     &   +                 RXRAT(   355 )      ! ACTP+NO3=ACO3+NO2+HCHO
     &   +    5.3800D-01 * RXRAT(   358 )      ! MACP+NO3=0.5380D+00*ACO3+...
     &   +    7.0000D-01 * RXRAT(   360 )      ! MVKP+NO3=0.7000D+00*ACO3+...
     &   +    3.5400D-01 * RXRAT(   394 )      ! ACRO=0.3540D+00*ACO3+...

c..Loss frequency of ACO3 ( not including ACO3 + ACO3 )
      L8 =                 RKI(   177 ) * YC ( NO2         )   ! ACO3+NO2=PAN
     &   +                 RKI(   204 ) * YC ( NO          )   ! ACO3+NO=MO2+NO2
     &   +                 RKI(   246 ) * YC ( HO2         )   ! ACO3+HO2=0.4400D+...
     &   +                 RKI(   283 ) * YC ( MO2         )   ! ACO3+MO2=0.9000D+...
     &   +                 RKI(   301 ) * YC ( ETHP        )   ! ACO3+ETHP=...
     &   +                 RKI(   302 ) * YC ( HC3P        )   ! ACO3+HC3P=...
     &   +                 RKI(   303 ) * YC ( HC5P        )   ! ACO3+HC5P=...
     &   +                 RKI(   304 ) * YC ( ETEP        )   ! ACO3+ETEP=...
     &   +                 RKI(   305 ) * YC ( OLTP        )   ! ACO3+OLTP=...
     &   +                 RKI(   306 ) * YC ( OLIP        )   ! ACO3+OLIP=...
     &   +                 RKI(   307 ) * YC ( BENP        )   ! ACO3+BENP=...
     &   +                 RKI(   308 ) * YC ( TOLP        )   ! ACO3+TOLP=...
     &   +                 RKI(   309 ) * YC ( XYMP        )   ! ACO3+XYMP=...
     &   +                 RKI(   310 ) * YC ( XYEP        )   ! ACO3+XYEP=...
     &   +                 RKI(   311 ) * YC ( ISOP        )   ! ACO3+ISOP=...
     &   +                 RKI(   312 ) * YC ( APIP1       )   ! ACO3+APIP1=...
     &   +                 RKI(   313 ) * YC ( APIP2       )   ! ACO3+APIP2=...
     &   +                 RKI(   314 ) * YC ( APINP1      )   ! ACO3+APINP1=...
     &   +                 RKI(   315 ) * YC ( APINP2      )   ! ACO3+APINP2=...
     &   +                 RKI(   316 ) * YC ( LIMP1       )   ! ACO3+LIMP1=...
     &   +                 RKI(   317 ) * YC ( LIMP2       )   ! ACO3+LIMP2=...
     &   +                 RKI(   318 ) * YC ( LIMNP1      )   ! ACO3+LIMNP1=...
     &   +                 RKI(   319 ) * YC ( LIMNP2      )   ! ACO3+LIMNP2=...
     &   +                 RKI(   321 ) * YC ( RCO3        )   ! ACO3+RCO3=MO2+ETHP
     &   +    5.0000D-01 * RKI(   322 ) * YC ( ACTP        )   ! ACO3+ACTP=...
     &   +                 RKI(   323 ) * YC ( MEKP        )   ! ACO3+MEKP=...
     &   +                 RKI(   324 ) * YC ( KETP        )   ! ACO3+KETP=...
     &   +    7.3100D-01 * RKI(   325 ) * YC ( MACP        )   ! ACO3+MACP=...
     &   +                 RKI(   326 ) * YC ( MCP         )   ! ACO3+MCP=NO2+...
     &   +                 RKI(   328 ) * YC ( UALP        )   ! ACO3+UALP=...
     &   +                 RKI(   329 ) * YC ( BALP        )   ! ACO3+BALP=MO2+BAL1
     &   +                 RKI(   330 ) * YC ( BAL1        )   ! ACO3+BAL1=MO2+BAL2
     &   +                 RKI(   331 ) * YC ( ADDC        )   ! ACO3+ADDC=...
     &   +                 RKI(   332 ) * YC ( MCTP        )   ! ACO3+MCTP=HO2+...
     &   +                 RKI(   333 ) * YC ( ORAP        )   ! ACO3+ORAP=MO2+GLY
     &   +                 RKI(   334 ) * YC ( OLNN        )   ! ACO3+OLNN=HO2+...
     &   +                 RKI(   335 ) * YC ( OLND        )   ! ACO3+OLND=...
     &   +                 RKI(   336 ) * YC ( ADCN        )   ! ACO3+ADCN=HO2+...
     &   +                 RKI(   337 ) * YC ( XO2         )   ! ACO3+XO2=MO2
     &   +                 RKI(   353 ) * YC ( NO3         )   ! ACO3+NO3=MO2+NO2
     &   +                 RKI(   400 ) * YC ( BDE13P      )   ! ACO3+BDE13P=...
     &   +                 RKI(   488 ) * YC ( VROCP6AROP  )   ! ACO3+VROCP6AROP=...
     &   +                 RKI(   494 ) * YC ( VROCP5AROP  )   ! ACO3+VROCP5AROP=...
     &   +                 RKI(   500 ) * YC ( NAPHP       )   ! ACO3+NAPHP=...

c..Loss frequency of PAN
      L9 =                 RKI(    37 )                        ! PAN=ACO3+NO2
     &   +                 RKI(    38 )                        ! PAN=MO2+NO3
     &   +                 RKI(   139 ) * YC ( HO          )   ! PAN+HO=XO2+NO3+HCHO
     &   +                 RKI(   178 )                        ! PAN=ACO3+NO2

c..K8_8, R8_9, and R9_8 terms
      K8_8  = RKI(   320 ) * DTC

      R8_9  = ( RKI(    37 )
     &      +   RKI(   178 ) ) * DTC 

      R9_8  = ( RKI(   177 ) * YC( NO2 ) ) * DTC 

c..Solution of quadratic equation to get ACO3 & PAN
      CMN = 1.0 + L9 * DTC
      A = 2.0D0 * K8_8 * CMN
      B = CMN * ( 1.0D0 + L8 * DTC ) - R8_9 * R9_8
      C = CMN * ( YC0( ACO3 ) + P8 * DTC ) +  R8_9 * YC0( PAN )

      Q = -0.5D0 * ( B + SIGN( 1.0D0, B ) * SQRT( B * B + 4.0D0 * A * C ) )

      YCP( ACO3 ) = MAX( Q / A , -C / Q  )

      YCP( PAN ) = ( YC0( PAN ) +  R9_8 * YCP( ACO3 ) ) / CMN

      RETURN

      END
